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We study the quantum phase transition from a super solid phase to a solid phase of p = 1/2 
for the extended Bose-Hubbard model on the honeycomb lattice using first principles Monte Carlo 
calculations. The motivation of our study is to quantitatively understand the impact of theoretical 
input, in particular the dynamical critical exponent z, in calculating the critical exponent v. Hence 
we have carried out four sets of simulations with ft = 2N 1/2 , f3 = SN 1/2 , f3 = N/2, and ft = N/4, 
respectively. Here j3 is the inverse temperature and N is the numbers of lattice sites used in the 
simulations. By applying data collapse to the observable superfluid density p„2 in the second spatial 
direction, we confirm that the transition is indeed governed by the superfluid-insulator universality 
class. However we find it is subtle to determine the precise location of the critical point. For 
example, while the critical chemical potential (p/V) c occurs at (p/V) c = 2.3239(3) for the data 
obtained using /3 = 2N 1 ^ 2 , the (p/V) c determined from the data simulated with /3 = N/2 is found 
to be (p/V) c = 2.3186(2). Further, while a good data collapse for p S 2N can be obtained with 
the data determined using (3 = N/4 in the simulations, a reasonable quality of data collapse for 
the same observable calculated from another set of simulations with /3 = 8N 1/I2 can hardly be 
reached. Surprisingly, assuming z for this phase transition is determined to be 2 first in a Monte 
Carlo calculation, then a high quality data collapse for p a iN can be achieved for (p/V) c ~ 2.3184 
and v ~ 0.7 using the data obtained with j3 — 8N 1 / 2 . Our results imply that one might need to 
reconsider the established phase diagrams of some models if the accurate location of the critical 
point is crucial in obtaining a conclusion. 

PACS numbers: 



I. INTRODUCTION 

Searching and investigating models in which one might observe a supersolid phase (SS) has been one of the central 
research interests in condensed matter physics recently. Indeed, it is reported numerically that several spin models 
and extended Bose-Hubbard models provide convincing evidence for the existence of a supersolid state, where the long 
range superfluid order and solid order coexist [H-[l2^. For example, a Monte Carlo investigation of bilayer spin-1/2 
Heisenberg model with an external uniform magnetic field on the square lattice demonstrates that a field-induced 
supersolid phase can be stabilized when the magnitude of the external magnetic field takes certain values [(| fll| . 
A supersolid phase is observed numerically by considering interacting bosons as well [lH5l l8Hl"ol. Il2| . Further, it is 
expected that SS and the quantum phase transitions out of SS might be realized experimentally using ultracold atoms 
on optical lattices [l3l - fT5| . Although rich phase diagrams have been obtained for both spin models and extended 
Bose-Hubbard models, a detailed investigation of the nature of the phase transitions out of SS is available only for 
bilayer spin-1/2 Heisenberg model on the square lattice. Further, despite the fact that the bilayer spin-1/2 Heisenberg 
model can be mapped into the hard-core Bose-Hubbard model, it is known that both the hard-core Bose-Hubbard 
model on the square and the honeycomb lattices do not exhibit a stable SS phase in the t/V-fi/V phase diagram at a 
fixed U/t [113. Here t is the nearest-neighbor hopping parameter, p, is the chemical potential, V is a nearest-neighbor 
^ \ repulsion, and U stands for an onsite repulsion. As a consequence, it will be interesting to carry out an investigation 
- - 1 on the nature of quantum phase transitions out of SS for the extended Bose-Hubbard model. Finally, considering the 
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many noticeable properties of graphene [17H20I ]. for which the underlying lattice is a honeycomb lattice, as well as 
the fact that quantum fluctuations are expected to be more relevant on the honeycomb lattice due to its coordination 
number, in this paper we study the phase transition from a SS phase to a solid state of p = 1/2 for the extended 
Bose-Hubbard model on the honeycomb lattice using first principles Monte Carlo simulations. Here p refers to the 
average number of bosons per lattice site. 



* fjjiang@ntnu.cdu.edu 
1 gcrbcru@itp.unibc.ch 



2 



The nature of quantum phase transitions out of a SS to a superfluid and a solid phases have been studied quanti- 
tatively using first principles unbiased quantum Monte Carlo method for the bilayer spin- 1/2 Hcisenberg model [TTj . 
Indeed it is demonstrated convincin gly that the critical theories determined from the Monte Carlo data are consis- 
tent with the theoretical prediction For example, by applying the technique of data collapse to the observable 
superfluid density p s , one reaches a perfect agreement between the numerical data and the theoretical prediction 
of superfluid- insulator universality class for the phase transition from a SS to a solid state [HI, E3 • To obtain the 
ground-state phase diagrams of the desired models using finite temperature Monte Carlo algorithms, in particular 
to quantitatively study the nature of the phase transitions in the phase diagrams, one useful strategy is to scale the 
inverse temperature (3 by the relation /? = cL z in the simulations. Here in addition to the dynamical critical exponent 
z which is already introduced earlier, L and c appearing in (3 = cL z are the box sizes used in the simulations and 
a constant, respectively. Surprisingly, despite the fact that the correct way to determine the critical theory for a 
second order phase transition is to employ the relation (3 = cL z in the simulations, several studies seem to simply use 
f3 = cL and ignore the impact of the dynamical critical exponent z when obtaining the ground-state phase diagrams. 
Whether the strategy of using (3 = cL instead of /3 = cL z in the Monte Carlo calculations has noticeable influence on 
determining the critical theory for a second order phase transition remains to explore. 

The motivation of our investigation on the quantum phase transition from a supcrsolid state to a solid state of 
p = 1/2 by varying p/V at fixed U/t = 20 and t/V = 0.16 in the parameter space for the extended Bosc-Hubbard 
model is twofold. First of all, we would like to determine the critical chemical potential (p/V) c as precise as possible 
since such a study is useful in calculating critical exponents such as {3 {13 /v) and r\. Secondly, since one should scale the 
inverse temperature (3 with the system size for the phase transition considered here due to the theoretical prediction 
z = 2, we would like to understand the impact of scaling j3 linearly with the linear length of the system on studying 
the corresponding critical theory of this transition. Indeed as we will demonstrate later, scaling (3 linearly with the 
system linear length either leads to poor data collapse for the obscrvables measured in this study, or one would arrive 
at a different critical theory than the expected one. 

This paper is organized as follows. After a brief description of the motivation behind this study, the extended Bosc- 
Hubbard model considered here as well as the observables measured in our Monte Carlo simulations are introduced. 
Follows that we present our numerical results. In particular, the corresponding critical point is determined with high 
precision by the method of data collapse. The subtlety of determining the critical theory for the transition from a SS 
phase to a solid state of p = 1/2 is demonstrated as well. Finally, a section is devoted to conclude our investigation. 

II. THE MODEL AND OBSERVABLES 

The extended Bose-Hubbard model considered in this study is given by 

H = -* E fa + b \ b ) + 1 E x ) 

+V J~] nj-rij - pi rij, (1) 

(id) 1 

where b\ (bi) are the bosonic creation (annihilation) operators at site i and rij is the occupation number at lattice site 
i. Further, U, V , t, and p appearing in Eq. ([T]) are defined as before. 

The honeycomb lattice with periodic spatial boundary conditions implemented in our simulations is depicted in 
figure 1. The dashed rectangle in figure 1, which contains 4 spins, is the elementary cell for building a periodic 
honeycomb lattice covering a rectangular area. For instance, the honeycomb lattice shown in figure 1 contains 3x3 
elementary cells. The lattice spacing a is the distance between two neighboring sites. The honeycomb lattice is not a 
Bravais lattice. Instead it consists of two triangular Bravais sub-lattices A and B (depicted by solid and open circles 
in figure 1). As a consequence, the momentum space of the honeycomb lattice is a doubly-covered Brillouin zone of 
the two triangular sub-lattices. 

Following earlier works in || [l(| , our primary interest for the model described by Eq. (TTJ) is to study the phase 
transitions from a supersolid phase to a solid phase of p = 1/2 and from a superfluid phase to a Mott insulator 
phase of p = 1 by varying the chemical potential p/V at fixed U/t = 20 and t/V = 0.16 in the parameter space. In 
particular wc would like to determine the critical chemical potentials (p/V) c as well as the correlation length critical 
exponents v and dynamical critical exponents z for these two phase transitions accurately To fulfill these tasks, 
we have measured the superfluid density p S i = (Wf)/f3 with i € {1,2} in our simulations and have focused on the 
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finite-size scaling of p S 2N z / 2 * . Here (W 2 ) is the winding number fluctuation in i-direction. 




FIG. 1: The periodic honeycomb lattice consisting of two triangular sub-lattices A and B, which are depicted by solid and 
open circles, respectively. The dashed rectangle is an elementary cell for building a periodic honeycomb lattice covering a 
rectangular area. 



III. DETERMINATION OF THE CRITICAL POINT 

To determine the location of the critical point in the parameter space fi/V, one useful technique is to study the 
finite-size scaling of certain obscrvables. For example, if the transition is second order, then near the transition, the 
observable p S iN z / 2 for i e {1, 2} should be described well by the following finite-size scaling ansatz 

N {j)=go{3N x ' 2v ,N*' 2 l[S\ (2) 

where On stands for p S iN z / 2 , j = (p c — p)/p c , and v is the critical exponent corresponding to the correlation length 
£. Finally go appearing above is a smooth function of the variables jN 1 / 2 " and N z l 2 /(3. In writing Eq. @, we have 
ignored explicitly the confluent correction to the scaling. Taking this correction into account, one concludes that for 
large enough (3 when the finite-temperature effects can be ignored, the curves of different N 1 / 2 for On, as functions of 
p/V, should have the tendency to intersect at critical point (p/V) c for large N 1 / 2 . Theoretically, it is predicted that 
the transition from a supersolid phase to a solid phase of p = 1/2 for this model is governed by the superfiuid-insulator 
universality class, namely one has v = 0.5, j3 = 0.5 and z = 2 [2l|, [22[. In the following we will apply the finite-size 
scaling formula, Eq. @, to the observable p s iL to determine (p,/V) c . In particular, we would like to examine whether 
the theoretical prediction is consistent with our Monte Carlo data. 



* Strictly speaking, the observable p s ; is defined by p a i = (Wf)/(pt). However since the parameter t is fixed to be t/V = 0.16 in our 
investigation, we will use p a i = {Wf)/f3 instead since the conclusions are not affected. 
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FIG. 2: Monte Carlo data of psiN as functions of p/V for the quantum phase transition from a supersolid state to a solid state 
of p = 1/2. For a given fixed N, the inverse temperature /3 used for this set of simulations is given by /3 = 2N 1 ^ 2 . 




FIG. 3: Data collapse for the observable p S 2N with 16 2 < N < 32 2 for the quantum phase transition from a supersolid state 
to a solid state of p — 1/2. The data is calculated using /3 = 2N 1 ^ 2 . In obtaining the figure, the critical exponent v is fixed to 
be the expected theoretical value v = 0.5. 



To study the critical theory for the phase transition of our central interest, we have performed large scale Monte 
Carlo simulations using the directed loop algorithms available in the ALPS library [23| . Without losing the generality, 
in our simulations we have fixed V to be 1.0 and have varied pi. To determine (p,/V) c using the finite-size scaling 
ansatz Eq. ([2]) , one needs to use large enough j3 so that the finite-temperature effects are negligible. For this purpose, 
we have carried out a trial simulation at pijV = 2.53 with N = 14 2 and we have found that one already obtains the 
zero-temperature result for p S 2 using j3 — 2 x 14 1 / 2 . Hence we use j3 = 27V 1 / 2 in other simulations as well. Notice the 
strategy of applying j3 = 2N 1 / 2 or similar ones for a fixed N 1 / 2 in the simulations was used in many studies exploring 
the phase diagrams of certain models in the literature. After determining the relation (3 = 2N 1 / 2 which allows one 
to access zero-temperature values for the observable p S 2, we have further carried out large scale simulations with N 
ranging from N = 10 2 to N = 32 2 . Figure [2] demonstrates the results of p S 2N as functions of p,/V for this set of 
simulations. The figure indicates that the phase transition is likely a second order transition because the curves of 
different N 1 / 2 have the tendency to intersect at a particular fi/V is the parameter space. Using fixed z = 2 and 
v = 0.5, the best result of data collapse for p S 2N is reached with (p/V) c = 2.3239(3) for the data of lattice sizes 
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16 2 < N < 32 2 (figure [3]). The quality of data collapse shown in figure [3] is not good, but acceptable. One might 
attribute the poor quality shown in figure [3] to a correction in Eq. ([5]) that is not taken into account in our analysis. 
At this stage, one would naturally conclude that our Monte Carlo data is consistent with the theoretical prediction, 
namely the critical exponent v and dynamical critical exponent z of the phase transition considered above is governed 
by v = 0.5 and z = 2. 
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FIG. 4: p s2 as a function of f3 at p/V = 2.41 and TV = 196. 



To make sure that we indeed obtain the ground-state properties of the model, we repeat above analysis by firstly 
determining the required (3 for reaching the ground-state value of p s2 at p/V = 2.41 and N = 14 2 . Surprisingly, we 
find that one has to use j3 ~ QN 1 / 2 in order to reach the zero-temperature value of p S 2 (figure 2]). Because of this 
observation, we have performed another set of simulations using j3 = 87V 1 / 2 . The results of p 8 2N as functions of /J./V 
for this new set of runs is shown in figurc[5] By comparing figures [2] and [5l one clearly observes a statistically difference 
between the critical chemical potentials calculated from these two set of data. To make the discrepancy between these 
two critical chemical potentials more transparent, we additionally simulating the model using p = 2N 1 / 2 in the range 
2.316 < fi/V < 2.32075 where (fj,/V) c for the second set of data (which are determined with f3 = 87V 1 / 2 ) is located. 
Figure [6] shows the results of p S 2-/V as functions of p/V for these new runs. No intersection between the curves of 
different TV 1 / 2 shown in figure [6] confirms our observation that the (p/V) c for these two set of data are statistically 
different. 
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FIG. 5: Monte Carlo data of p a2 N as functions of p/V for the quantum phase transition from a supersolid state to a solid state 
of p — 1/2. For a given fixed N, the inverse temperature /3 used for this set of simulations is given by (3 = 8N 1 ' 2 . 
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FIG. 6: Monte Carlo data of p a2 N as functions of p,/V in the range 2.316 < ft < 2.32075 with 10 2 < N < 32 2 . For a given 
fixed N , the inverse temperature f3 used for this set of simulations is given by /3 = 2N 1 ^ 2 . 



After demonstrating that the critical chemical potentials determined from the sets of data obtained using /? = 2A 1//2 
and fi = 8A 1 / 2 are statistically different, let us return to the analysis of the second set of p S 2 data. Yet another 
surprise we find is that a good data collapse for p S 2N can hardly be achieved if data points of small A are included in 
the analysis. Only with data of large A, namely A > 24 2 , a reasonable data collapse for p s2 A with a fixed v = 0.5 can 
be obtained given that (p/V) c = 2.3185(2) (figure[7j)- The results presented so far in this study raise an interesting 
question, namely what values of (3 should be used in our simulations in order to obtain the ground-state values for the 
observable p S 2- In other words, should one perform another set of simulations with even larger values of /?? Actually 
by examining the relevant finite-size scaling ansatz Eq. ([2|) carefully, one would realize that the correct strategy is to 
scale f3 with N. Figure [8] shows the results of p S 2N obtained with j3 = N/2 as functions of p,/V. The tendency of 
crossing between curves of different A 1 / 2 shown in figure [5] is much stronger compared to those presented in figures 
[2]and[5l Finally the quality of data collapse for p S 2N demonstrated in figure [9] which are obtained from the set of 
simulations using /3 = A/2 is also much better than those found in figures [3] and [7] Notice the (p,/V) c determined 
from the data simulated with (3 = A/2 is given by (p,/V) c = 2.3186(2) which agrees with (p,/V) c = 2.3185(2), but is 
statistically different from (p,/V) c = 2.3239(3) obtained earlier from the data determined using /3 = 2A 1 / 2 . 
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FIG. 7: Data collapse of the observable p s %N with 24 2 < N < 32 2 for the quantum phase transition from a supersolid state to 
a solid state of p = 1/2. The data is calculated using 13 = 8N 1/2 . In obtaining the fi gure, the critical exponent v is fixed to be 
the expected theoretical value v — 0.5. 




FIG. 8: Monte Carlo data of p S 2N as functions of fi/V for the quantum phase transition from a supersolid state to a solid state 
of p — 1/2. For a given fixed N, the inverse temperature /3 used for this set of simulations is given by /3 = N/2. 
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FIG. 9: Data collapse of the observable p S 2N with 16 2 < N < 32 2 for the quantum phase transition from a supersolid state to 
a solid state of p = 1/2. The data is calculated using /3 = N/2. In obtaining the figure, the critical exponent v is fixed to be 
the expected theoretical value v — 0.5. 



Of course, one might argue that since N/2 > 8N 1 / 2 in the range of N used in the simulations and analysis, the 
better results shown in figures [8] and [9J than those in figures El El El and [7j is simply because one reaches the ground- 
state values of p S 2 for the set of data determined using B = N/2 and the p S 2 data points calculated with 8 = 8N 1 / 2 
still receive finite-tcmpcraturc effects. To rule out such a possibility, we have performed another set of simulations 
using 8 = N/4 with which the condition N/4 < 8N 1 / 2 is satisfied for N 1 / 2 = 18,20, ...,32. Figures HU1 and fTTI show 
the results of p S 2N with 18 2 < N < 32 2 for these new simulations using B = N/4 as functions of p/V and the 
corresponding results of data collapse. In obtaining figure fTTI the critical exponent v is fixed to v = 0.5 as before. 
Again the quality of crossing and data collapse seen in figures [JO] and QTJ are much better than those in figures O and 
[71 To make the comparison on the same footing, figure [T2] shows the results of data collapse for p S 2-/V determined 
using B = 87V 1 / 2 with 18 2 < N < 32 2 . Fi gures [TTJ and [JJ] clearly indicate that the quality of data collapse in figure fTTI 
is much better than that of figure [TSJ Interestingly, assuming z = 2 is determined first in a Monte Carlo calculation 
before one performs the large scale simulations for obtaining p S 2 using B — 87V 1 / 2 , then a good data collapse can 
be reached with (p,/V) c = 2.3184 and v = 0.7 using the data of p s2 N(N > 24 2 ) calculated with B = 87V 1 / 2 (figure 
[TBI) . While {p/V) c = 2.3184 is consistent with (p/V) c = 2.3185(2) and {p,/V) c = 2.3186(2), the v = 0.7 we find is 
significantly different from the theoretical expectation v = 0.5. Without the theoretical input v = 0.5, one might be 
misled by figure [T3] to conclude an unconventional phase transition is observed for this model. Our results presented 
in this study clearly imply the importance of using the correct relation for B and N 1 / 2 in investigating the critical 
theory of a second order phase transition. 
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FIG. 10: Monte Carlo data of psiN as functions of p/V for the quantum phase transition from a supersolid state to a solid 
state of p — 1/2. For a given fixed N, the inverse temperature /3 used for this set of simulations is given by /3 = N/4. 
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FIG. 11: Data collapse of the observable p S 2N with 18 2 < N < 32 2 for the quantum phase transition from a supersolid state 
to a solid state of p — 1/2. The data is calculated using /3 = N/4. In obtaining the figure, the critical exponent v is fixed to 
be the expected theoretical value v — 0.5. 
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FIG. 12: Data collapse of the observable p S 2N with 18 2 < N < 32 2 for the quantum phase transition from a supersolid state 
to a solid state of p = 1/2. The data is calculated using f] = 8iV 1/2 . In obtaining the fi gure, the critical exponent v is fixed to 
be the expected theoretical value v — 0.5. 



2.2 r 

2 j- 
1.8 f- 

1.6r 

1.4 =- 

« ih 

°- 0.8 H 

0.6- 
0.4- 
0.2- 
0^ 
-0.2^ 
-0.4^ 



0i/V) c = 2.3184 

N >= 576 
v = 0.7 



lJ_i 



J_l 



lJ_i 



-0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0.05 0.1 0.15 0.2 



04i c -1)N 



1/2v 



FIG. 13: Data collapse of the observable p s iN with 24 2 < N < 32 2 for the quantum phase transition from a supersolid state 
to a solid state of p = 1/2. The data is calculated using ft = 87V 1//2 . In obtaining the figure, the critical chemical potential p c 
and the critical exponent v are fixed to be p, c = 2.3184 and v = 0.7, respectively. 



In additional to the quantum phase transition from a supersolid phase to a solid phase of p = 1/2 for the extended 
Bose-Hubbard model, we have studied the critical theory of the phase transition from a superfluid phase to a Mott 
insulator phase of p = 1 as well. Inspired by the subtlety we observed earlier, we use (3 = N/2.5 for this new 
investigation since the dynamical critical exponent z for this transition is predicted to be 2 theoretically. Indeed 
using v = 0.5 which is the expected theoretical value for v, a good data collapse is reached for p S 2N provided that 
{p/V) c = 3.6378(2) and N > 16 2 (figure [14]). In other words, our Monte Carlo data is fully compatible with the 
predicted universality class for the transition from a superfluid phase to a Mott insulator phase of p = 1 . 
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FIG. 14: Data collapse of the observable p B 2N with 18 2 < N < 32 2 for the quantum phase transition from a superfluid state 
to a Mott insulator state of p = 1. The data is calculated using ft = N/2.5. In obtaining the figure, the critical exponent v is 
fixed to be the expected theoretical value v = 0.5. 



IV. DISCUSSIONS AND CONCLUSIONS 



In this paper, we use first principles Monte Carlo methods to study the quantum phase transitions from a supcrsolid 
phase to a solid phase of p — 1/2 and from a superfluid phase to a Mott insulator phase of p = 1 by varying the 
chemical potential p,/V at fixed U/t = 20 and t/V = 0.16 in the parameter space for the extended Bosc-Hubbard 
model. We confirm that our Monte Carlo results for both phase transitions arc fully compatible with the corresponding 
theoretical prediction. Specifically, we obtain a good quality of data collapse for p S 2N with fixed v — 0.5 and z = 2 
in the analysis. On the other hand, we find it is subtle to determine the correct critical chemical potential (n/V) c for 
the phase transition from a supersolid phase to a soild phase of p = 1/2. For example, while the (p/V) c calculated 
from a set of p S 2 data obtained using ft = 2N 1 / 2 is given by {jj,/V) c = 2.3289(2), another set of p S 2 data determined 
from the simulations with ft = N/2 leads to (p/V) c = 2.3186(2). The deviation of these two values of (p/V) c found 
from different strategies of scaling the inverse temperature ft with the system sizes is statistically significant. The 
later result for (p/V) c should be more reliable considering the fact that the dynamical critical exponent z is 2. We 
also demonstrate that using the incorrect scaling of ft, namely one scales ft linearly with the system linear length 
N 1 / 2 , would either lead to bad quality of data collapse for the observables considered here, or one would arrive at 
a different critical theory for the phase transition investigated in this study. Despite the fact that the principles 
behind the technique of data collapse is exact, it might mislead one to conclude an observation of an unconventional 
phase transition. Our results of (p/V) c — 2.3184 and v = 0.7 obtained by carrying out the data collapse for p S 2-/V 
determined from the simulations with ft = 8./V 1 / 2 is a good example of such a scenario. When applying the technique 
of data collapse to study the critical theory of a second order phase transition, it might be useful to employ this 
technique in a more rigorous manner. At this stage, it is a surprise that the motivation behind our investigation, 
namely whether scaling ft linearly with the system linear length for studying a second order phase transition with 
z — 2 is an appropriate strategy, is ignored in some studies of phase diagrams for certain models in the literature. 
To re-examine the corresponding critical points seems to be the required step in order to confirm the quantitative 
correctness of the established phase diagrams of some spin and bose models. Our results does not necessarily indicate 
that the phase diagram of the bilaycr spin-1/2 Heisenberg model in an external magnetic field studied in [TT| and 
the phase diagram of extended Bose- Hubbard model mapped out in [Tfj| are not quantitatively correct. However 
one might need to reconsider the phase diagrams of these models if the precise location of a critical point in the 
parameter space is crucial in drawing a conclusion. Finally, based on the results presented in this paper, in particular 
the unexpected u = 0.7 obtained from the data p S 2N calculated with ft = 8N 1 / 2 , it will be important to determine 
the dynamical critical exponent z first when investigating the critical theory for a second order phase transition which 
might be described by a new universality class and is poor known theoretically. 
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